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1.  FOREWORD 


One  of  important  goals  of  the  demining  mission  is  a  reliable  differentiation  between  land 
mines  and  clutter,  which  would  decrease  the  false  alarm  rate.  This  project  can  be  considered 
as  a  new  approach  to  this  issue,  via  a  new  method  for  solutions  of  inverse  problems.  It  is 
natural,  therefore,  to  test  such  an  approach  on  a  mathematical  model  first  in  order  to  see 
its  main  features.  While  being  at  an  Army’s  demonstration  of  a  demining  process  at  AP 
Hill  (VA)  in  1998,  we  observed  that  the  process  of  detection  and  identification  of  land  mines 
by  a  hand-held  GPR  can  be  subdivided  into  two  subprocesses,  which,  using  an  analog  with 
medical  imaging,  we  call:  (i)  screening  and  (ii)  diagnostics.  In  the  screening  procedure 
one  identifies  horizontal  coordinates  of  possible  mine-like  targets,  many  of  which  might  be 
clutter.  The  screening  procedure  should  be  very  rapid.  The  diagnostic  procedure,  however, 
is  slower.  This  is  a  procedure  of  differentiating  between  mines  and  clutter.  Another  term 
for  this  is  “classification”  procedure.  All  targets  being  detected  on  the  first  stage  should 
be  checked  out  again.  Finally,  all  targets,  being  identified  as  possible  land  mines  should 
be  removed  from  soil  by  soldiers.  Since  the  diagnostic  procedure  usually  takes  many  hours, 
even  for  a  small  area,  it  would  be  quite  helpful  to  develop  new  methods,  which  would  speed 
it  up  very  essentially,  while  still  maintaining  high  reliability. 

Prom  this  prospective,  the  goal  of  this  project  was  to  develop  a  new  mathematical  ap¬ 
proach  of  an  inverse  problem,  which  potentially  might  lead  to  a  successful  diagnostic  proce¬ 
dure  for  land  mines.  The  input  data  for  an  inverse  problems  are  measurements  by  a  GPR  of 
a  back  reflected  electric  signal  on  an  array  of  antennas,  using  a  frequency  sweep.  The  output 
of  the  solution  of  the  inverse  problem  is  the  vector  consisting  of  coordinates  of  a  target(s), 
as  well  as,  dielectric  permittivity  e  and  electric  conductivity  a  within  a  target(s).  On  the 
next  stage,  this  vector  is  supposed  to  be  used  as  an  input  for  the  classification  procedure. 
Since  we  are  interested  in  imaging  of  plastic  antipersonnel  mines  with  a  low  depth  of  burial 
ranging  from  1  cm  to  10  cm,  we  work  with  a  high  frequencies  penetrating  up  to  this  depth. 
The  range  of  these  frequencies  is  from  0.5  GHz  to  3  GHz. 

The  solution  of  the  corresponding  inverse  problem,  however,  is  quite  a  challenging  issue 
in  its  own  right.  This  issue  was  addressed  in  this  project.  Specifically,  a  new  method  for  the 
corresponding  inverse  problem  was  developed,  which  is  a  second  generation  of  the  so-called 
Elliptic  Systems  Method  (ESM),  being  previously  developed  by  the  Co-PIs  for  Diffusion 
Tomography  with  medical  applications  [10,  12-14].  (In  this  report  references  in  bold  face 
letters  refer  to  Bibliography  section.  References  like  [1]  refer  to  section  fitted  ’’Listing  of  all 
publications  and  technical  reports  supported  under  this  grant.”)  The  new  method  provides 
accurate  estimates  of  the  above  parameters  in  about  six  minutes  of  CPU  time  on  SGI  Origin 
200  with  one  processor.  This  time  frame  is  realistic  for  the  diagnostic  purpose.  It  can  also 
be  decreased  if  using  several  processors.  The  key  reason  for  such  a  small  time  frame  is  that 
in  this  new  approach  the  resulting  linear  system  has  a  differential,  rather  than  conventional 
integral  form.  This  means,  in  turn,  that  the  matrices  to  be  inverted  are  sparse,  rather  than 
full,  the  very  property,  which  enables  one  to  apply  modern  techniques  of  Numerical  Linear 
Algebra  for  their  inversion.  We  want  to  point  out  that  almost  all  conventional  methods  of 


1 


solutions  of  inverse  problems  are  relying  on  solutions  of  ill-posed  integral  equations,  which 
inevitably  lead  to  large  full  matrices  to  be  inverted,  thus  requiring  a  significant  CPU  time. 

Since  methods  of  integral  equations  are  quite  time  consuming,  as  compared  with  ours, 
we  compare  our  technique  with  the  best  competing  approach,  being  initially  developed  by  a 
well  known  German  expert  in  inverse  problem  F.  Natterer  six  years  ago  [19];  also  see  follow 
up  publications  [5,6].  Thus,  we  modified  our  algorithm  accordingly.  Naturally,  the  modified 
variant  is  also  based  on  resulting  differential,  rather  than  integral  operators.  In  the  view  of 
possible  applications  to  the  experimental  data  an  important  advantage  of  the  second  version 
is  that  the  differentiation  of  the  data  with  respect  to  the  frequency  is  not  required,  unlike  the 
first  version.  On  the  other  hand,  while  both  versions  have  a  similar  performance  in  terms 
of  timing  and  locations  of  targets,  the  first  one  still  often  provides  better  values  of  electrical 
parameters  within  targets,  which  is  the  key  for  the  diagnostic  goal.  Since  the  development  of 
computational  tols  of  the  project  was  very  time  consuming,  we  did  not  have  sufficient  time 
left  to  work  on  the  experimental  data.  Our  modified  algorithm  is  an  advanced  version  of 
Natterer’s  method,  since  we  calculate  the  resulting  matrices  precisely,  whereas  only  diagonal 
elements  of  matrices  are  counted  in  [5,6,19].  Thus,  if  using  only  those  elements  (i.e.,  if 
literally  following  the  idea  of  Natterer),  images  by  the  second  method  would  be  much  worse 
than  by  the  first  one. 

Although  the  topic  of  the  project  is  inverse  problems,  one  must  always  possess  a  rapid 
algorithm  for  solutions  of  corresponding  forward  problems,  when  working  with  the  inverse. 
However,  when  starting  the  work  on  this  project,  we  discovered  that  such  an  algorithm 
simply  did  not  exist  at  that  time,  for  the  high  frequency  range  we  used.  Therefore,  the  very 
first,  although  axillary  task  was  to  develop  such  a  method. 

Thus,  major  tasks  of  this  project  are  ones  listed  below.  Work  on  each  of  these  tasks 
required  approximately  one  year. 


Task  1  (axillary) 

Development  of  an  efficient  numerical  method  for  the  solution  of  the  forward  GPR  prob¬ 
lem  for  the  high  frequencies. 


Task  2 

Development  of  a  second  generation  of  the  ESM  for  the  solution  of  the  GPR  inverse 
problem. 


Task  3 

Development  of  the  modified  version  of  the  method  of  Task  2  and  comparison. 

All  these  three  tasks  were  successfully  completed.  Summary  of  results  in  given  below. 
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4.  STATEMENT  OF  THE  PROBLEM  UNDER  STUDY 

4.1.  Partial  Differential  Equation 

The  GPR  signal  is  modeled  here  as  a  polarized  electric  plane  wave  Eq  =  (0,0,  exp 
{iu.y/fIeo  •  0:2))  •  exp{—iu}t)  propagating  along  the  positive  direction  of  a;2-axis  in  the  half 
space  {x2  <  0}.  Here  u  =  27r/  is  the  angular  frequency  of  the  signal,  //  =  47r  •  10“^  Henry/m 
is  the  magnetic  permeability  and  Eq  =  8.854  •  10"^^  Farad/m  is  the  dielectric  permittivity  of 
free  space.  It  is  assumed  that  {x2  <  0}  is  air  and  {x2  >  0}  is  ground,  where  mine-like  tar¬ 
gets  are  located.  All  functions  below  depend  only  on  two  spatial  variables  {x\,X2)  =  x.  Let 
E{x,(jj)  =  (0, 0,«(a:,a;))  •  exp{—i(jjt)  be  the  electric  field.  Then  the  following  Helmholtz-like 
PDE  for  the  function  u(x,  u)  can  be  derived  from  Maxwell’s  system  [3] 

+  k^{x,u)u  =  0.  (4.1) 

Here  the  function  has  the  form  k'^{x,uj)  =  ui^iJ,e{x)  +  iojfia{x),  where  e(x)  and  a{x) 

are  respectively  the  dielectric  permittivity  and  the  electric  conductivity  of  the  medium,  and 
/X  is  the  magnetic  permeability,  which  is  assumed  to  be  constant  everywhere  and  equal  to  its 
value  in  firee  space,  /x  =  const.  We  also  assume  that  e  =  Sq  in  air. 

4.2.  Ranges  of  Parameters 

It  is  useful  first  to  establish  ranges  of  parameters  in  the  PDE  (4.1).  All  units  below  are 

given  in  SI  system.  The  frequency  of  the  signal  /  =  —  is  between  0.5  GHz  and  3  GHz,  i.e., 

27r 

/  6  (0.5,3)  •  10® — .  Let  e  =  £r£o  where  Sr  is  the  relative  dielectric  constant.  In  air  Sr  =  1 
sec 

and  cr  =  0.  We  introduce  the  so-called  ’’loss  tangent”  as  [9] 

tan(5)  =  (4.2) 

(jJE 

Then 

k^  =  +  i  •  tan(5)].  (4.3) 

d 

We  assume  that  the  loss  tangent  does  not  depend  on  u;,  i.e.,  —  [tan(5)]  =  0.  This  condition 

ULO 

is  a  requirement  for  the  presented  imaging  algorithm.  It  is  satisfied  with  sufficient  accuracy 
in  many  practical  scenarios  of  land  mine  detection. 

The  approximate  values  of  the  parameters  tan((5),  A;  and  the  wavelength  A  = 
27r/Re(A:)  for  different  soil  moistures  as  well  as  for  trinitrotoluene  (TNT)  are  given  in  Table 
1  for  the  frequency  /  =  1  GHz.  In  this  table  we  use  the  data  of  [7] 

Table  1.  Approximte  Values  of  er,tan(<5),  A:^,  A:  and  A  for 


Different  Soil  Moistures  and  TNT  at  /  =  1  GHz 


Medium 

B 

tan((5) 

A:2 

Qj^HHIIIIIU 

A  [cm] 

Air 

1 

0 

439 

.2 

21 

30 

Dry  soil 

30 

0.025 

1273  +  ^  31 

35.7 -1-z- 0.43 

17 

Wet  soil,  5%  moisture 

4 

0.22 

1756  -H  ?  •  395 

42 +  i-  4.7 

15 

TNT 

2.86 

0.0018 

1256  -1-  i  •  2.26 

35.4  -h  i  •  0.03 

17.7 
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4.3.  Statement  of  the  Forward  Problem 

We  assume  that  the  electrical  parameters  e  and  a  have  constant  background  values 
everywhere  in  the  ground,  except  in  the  mine-like  targets,  whose  sizes  are  small,  as  compared 
with  the  size  of  the  region  of  interest.  Let  ko  =  ko{x2,u>)  be  the  function  k  in  (4.3)  for  the 
background  medium.  Then  this  function  has  a  discontinuity  on  the  air/ground  interface. 


f  u^fieo,  for  X2  <0 

\  u)^fieoSr[l  +  i  •  tan(^)],  for  X2  >  0. 


(4.4a) 


Let  uo  =  tto(x2,a;)  be  the  solution  of  the  PDE  (4.1),  which  corresponds  to  the  initial  plane 
wave  without  targets  present.  Then  hq  consists  of  the  initial,  reflected,  and  transmitted 
plane  waves  [9] , 

^ikoX2+Riko)e-*>’o^2  ^  <  0 

T(A:o)e*'=“^^  for  Xz  >  0, 

where  R{ko)  and  T(A:o)  are  the  reflection  and  transmission  coefficients  given  by 


Uo 


=  { 


(4.4b) 


R{ko)  = 


kp  ~  k^ 
kp  "b 


Tikp)  =  ^ 


2/uo 


kp  +k^' 


(4.4c) 


Here  k^  and  kp  are  the  values  of  kp  for  0:2  <  0  and  0:2  >  0  respectively.  The  presence  of 
these  coefficients  ensures  the  continuity  of  the  function  up  together  with  its  first  derivatives 
at  {x2  =  0}. 

We  seek  a  solution  of  the  equation  (4.1)  in  the  form  u  =  up-\-v,  where  the  function  v  {x,  oj) 
represents  the  wave  scattered  by  mine-like  targets  with  compact  supports  in  =  {x2  >  0}. 
Hence,  this  function  satisfies  the  following  PDE 


+  k^v  =  —g,  x  G  R^, 


(4.5) 


where  outside  of  the  targets  and 

g{x,oj)  =  I  (jj 


outside  targets 
(A:^  —  fco)'Uo,  inside  targets. 


In  addition,  we  impose  Sommerfeld  radiation  boundary  conditions  at  infinity 

]im^/r  (^  -  ikpv]  =  0, 

r-KX  \or  ) 


(4.6) 


where  r  =  \Jx\  -f  a;^,  Im(A:o)  >  0  and  the  limit  holds  uniformly  in  all  directions.  The 
uniqueness  and  existence  of  the  solution  of  the  problem  (2.6),  (2.7)  was  proven  in  [3]  for 
the  3-D  case.  It  follows  from  the  proof  in  [3]  that  similar  results  hold  in  the  2-D  case 
for  functions  v  G  with  s  >  1/2,  where  =  {cu  :  G  |q:|  <  1}  and 

IJ2,s^r2)  =  ;  (1  +  \x2\^y^'^v  G  £2(7?^)}.  Likewise,  if  D  is  any  bounded  domain  such  that 

either  D  C  R\  or  D  C  R!L,  then  v  G  H‘^{D).  It  follows  also  from  [3]  that  if  in  (4.4a) 
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tan((J)  >  0,  then  for  0:2  >  0  the  functions  v{x,ijj),Vx^{x,uj)  and  Vx2{x,(jj)  decay  exponentially 
as  a;  — »•  00.  The  function  {k'^—kQ){x)  is  assumed  to  be  bounded  and  to  have  compact  support 
in  {x2  >  7}  for  some  positive  constant  7. 

It  is  natural  to  consider  in  practical  computations  a  bounded  domain  Gi,  which  is  ob¬ 
tained  by  a  cut-off  of  the  infinite  space  R^.  For  the  numerical  solution  of  the  forward  problem, 
we  consider  a  square  Gi  =  \x<i\  <  L).  In  this  case  the  condition  (4.6)  is  replaced  with 

Vxi  T  ihv\xi=±L  =  0,  Vx^  T  ikov\x2=±L  =  0.  (4.7) 

So,  below  we  will  always  assume  that  the  boundary  value  problem  (4.5),  (4.7)  has  an  unique 
solution  V  e  H^{Gl)-  By  the  well  known  results  for  elliptic  equations,  this  implies,  in  turn, 
that  V  e  H‘^{Gl  n  {x2  >  0})  and  v  €  H^{Gl  H  {x2  <  0}).  A  natural  question  to  ask  would 
be  about  the  influence  of  the  value  of  the  cut-off  constant  L  on  the  resulting  solution  v.  It 
was  shown  numerically  in  [1]  that  if  the  targets  are  located  ’’well  within”  the  square  Gl  (i.e., 
far  from  the  boundaries),  then,  for  the  range  of  parameters  listed  in  Table  1,  the  resulting 
value  of  the  function  v{x,uj)  for  points  x  located  near  the  air/ground  interface  {x2  =  0} 
is  independent  of  L  as  long  as  L  >  53  cm.  For  this  reason,  for  the  solution  of  the  inverse 
problem  we  chose  a  smaller  rectangle  C  Gl,Q,  C  /?+.  First  we  generate  the  data  f|i2=o 
and  Vx2 \x2=o  for  the  inverse  problem  using  the  solution  of  the  forward  problem  in  the  domain 
Gl  with  L  =  1.5  m.  To  solve  the  inverse  problem,  we  use  the  rectangle 

=  {x  =  (xi, X2)  :  |xi|  <  Li  =  0.6  m,  0  <  X2  <  T2  =  0.4  m}. 

In  doing  this,  we  use  the  following  boundary  conditions  for  the  function  v  on  the  side  and 
top  boundaries  of  : 

Vxi  T  ikv\xi=±Li  ~  0  (4.8a) 

Vx2  -  ikv\x2=L2  =  0.  (4.8b) 

As  to  the  boundary  conditions  on  the  bottom  boundary  v\x2^o  ^'^id  Vx2\x2^ot  they  are  used 
as  an  input  data  for  the  inverse  problem  and,  thus  are  taken  from  the  solution  of  the  forward 
problem  in  the  domain  Gl- 

4.4.  Statement  of  the  Inverse  Problem 

Let  £1  and  tan(5i)  be  the  values  of  the  parameters  e  and  tan(5)  everywhere  in  the  ground, 
except  for  the  mine-like  targets.  Then 

e(x)  =  £1  +  /ie(x),  tan(5)  =  tan((5i)  -I-  ha{x),  (4.9) 

where  the  perturbations  h^{x)  and  ha{x)  are  due  to  the  presence  of  the  mine-like  targets. 
Hence,  the  determination  of  these  functions  would  yield  both  the  locations  of  these  targets 
and  the  values  of  the  electrical  parameters  within  them.  Since  the  loss  tangent  does  not 
depend  on  the  frequency  u),  we  introduce  a  perturbation  h{x)  of  the  background  coefficient 
kl(x  ,uj)  as  a  ’’whole,” 

^  k^{x,uj) -kl{x,uj) 

— 


(4.10) 


7 


Hence, 


k‘^{x,u))  =  A:o[l  +  h{r)]  =  +  h{x)],  (4.11) 

where  =  kl/uj^.  By  (4.9)-(4.11)  the  function  h{x)  can  be  obtained  from  the  functions 
he{x)  and  ha{x)  through  the  following  transformation: 

_  he  +  ijeihg  +  he  tan((^i)  +  hehg] 
ei[l  +  itan(5i)] 

Once  the  complex  valued  function  h{x)  has  been  obtained  from  the  solution  of  the  inverse 
problem,  one  can  recover  the  perturbations  of  the  physical  parameters  by  using  the  formulas: 

he  =  ei[Re(/i)  —  tan(5i)  •  Im(/i)], 

Im(h)  •  [1  +  tan^(5i)] 

1  +  Re(h)  —  tan((5i)  •  lm{h) 

Inverse  Problem.  Given  functions  ip{xi,u)),'^{xi,w)  defined  as 

(p{xi,oj)  =  v\x2=o,  (4.12a) 

ip{xi,uj)  =  Vx2\x2=o,  (4.12b) 

for  Xi  e  (— Li,Li),  u  G  (wniin)<^max))  determine  the  perturbation  function  h{x). 

Here  (wmim^'^max)  is  the  available  frequency  band,  over  which  measurements  are  per¬ 
formed.  We  assume  that  measurements  are  performed  at  points  on  a  certain  interval 
(— Li,Li)  of  the  line  {x^  =  0}  located  on  the  air/ground  interface.  To  explain  a  possi¬ 
ble  way  to  evaluate  the  normal  derivative  in  (4.12b),  we  observe  that  if  the  function 
is  given,  one  can  uniquely  solve  the  boundary  value  problem  (4.5),  (4.8),  (4.12a)  in  the  air, 
i.e.,  for  {x2  <  0}  D  {|xi|  <  Li},  where  no  targets  are  present.  This  is  sufficient  to  determine 
the  function  V’-  Another  option  to  obtain  the  normal  derivative  Vx2  would  be  to  measure 
not  only  the  third  component  of  the  electric  field  E,  but  the  first  component  Hi  of  the 
magnetic  field  H  as  well,  since  Maxwell’s  system  implies  in  our  case  E3X2  =  iojpHi. 

The  resulting  PDE  for  the  function  v  is 

+  ufiklgV  4-  J^kl^h{x)v  =  —ufiklgh{x)uQ,  (4-13) 

where  h{x)  is  a  bounded  function  with  compact  support  in  and  the  function  uq  is  given  by 
(4.4b)  (for  X2  >  0).  We  will  assume  that  the  medium  of  interest  is  basically  homogeneous, 
except  for  a  few  mine-like  targets  whose  sizes  are  small  compared  to  the  size  of  fi.  This 
suggests  the  related  assumption  that  ^  Ikllizfn)’  function  s{x)  =  1. 

Note  that  function  v  depends  nonlinearly  on  the  function  h.  However,  we  will  consider  a 
linearized  inverse  problem,  assuming  that  perturbations  are  small,  as  compared  with  the 
background.  Hence,  linearization  leads  to  dropping  the  h{x)v  term  in  (4.13).  This  approach 
was  used  previously  in  publications  about  the  ESM  [10,  12-14].  The  assumption  about  the 
linearization  can  actually  be  relaxed  if  using  Newton-like  updates  [10,  12-14].  However, 
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limited  testing  of  those  updates  for  the  range  of  parameters  we  used  did  not  show  significant 
improvement  of  the  images. 

Thus,  below  we  will  consider  the  inverse  problem  only  for  the  linearized,  with  respect  to 
h,  equation  (4.13),  while  the  data  simulation  for  the  forward  problems  will  still  be  done  for 
the  original  equation  (4.13).  Linearization  of  (4.13)  with  respect  to  h  leads  to 

+  uPkQgV  =  —u?klJi{x)uo  (4-14) 

5.  SUMMARY  OF  THE  MOST  IMPORTANT  RESULTS. 

In  this  section  we  discuss  only  results  of  tasks  1-3,  which  are  the  most  important  ones 
in  this  project.  These  results  were  obtained  by  M.  V.  Klibanov,  Yu.  A.  Gryazin  and  T.  R. 
Lucas.  It  should  be  noted  that  some  results  of  a  lesser  importance  were  also  obtained  by 
other  members  of  the  scientific  personnel :  T.  P.  Weldon,  V.  J.  Patel  [4,6]  and  S.  Pamyatnikh. 
However,  these  results  are  not  discussed  here. 

5.1.  TASK  1  (AXILLARY). 

DEVELOPMENT  OF  AN  EFFICIENT  NUMERICAL  METHOD  FOR  THE 
SOLUTION  OF  THE  FORWARD  GPR  PROBLEM  FOR  HIGH 

FREQUENCIES 

Since  this  is  axillary  (although  an  important)  task,  we  will  provide  a  rather  brief  descrip¬ 
tion  of  this  development  referring  to  [1,3]  for  details. 

As  it  was  pointed  out  above,  the  very  first  difficulty  we  faced  in  this  project  was  the 
absence  of  a  rapid  algorithm  for  the  solution  of  the  forward  problem  for  small  wavelengths 
A  we  used.  For  example,  as  it  follows  from  Table  1,  for  the  frequencies  /  €  (0.5,3)  GHz  the 

.  J  (5, 35)  cm,  in  the  ground 
^  (  (10, 60)  cm,  in  the  air 

This  range  of  wavelengths  significantly  affects  the  grid  size  in  the  Finite  Differences  Solution 
for  both  the  forward  and  the  inverse  problem.  In  order  to  calculate  the  forward  problem 
accurately,  one  should  use  at  least  10  grid  points  per  wavelength.  Suppose,  for  example  that 
one  wants  to  calculate  the  function  u  in  a  square  region  of  2  m  x  2  m.  Then,  because  of  (5.1), 
this  would  mean  that  one  should  use  at  least  400  x  400  grid  for  A  =  5  cm.  When  analyzing 
this  problem,  we  quickly  realized  that  standard  direct  solve  routines,  such  as  LAPACK,  for 
example,  are  inapplicable  here,  since  they  provide  data  generation  for  the  inverse  problem 
in  about  6  hours  of  CPU  time  on  Silicon  Graphics  Origin  200  (SGI).  This  was  a  strong 
motivation  for  us  to  develop  a  new  efficient  algorithm  for  the  forward  problem. 

The  idea  is  to  develop  a  high  quality  Krylov  subspace  based  method  (GMRES).  First, 
we  approximate  the  differential  operator  of  the  Helmholtz  equation  (4.5)  in  the  square  Gl 
with  finite  differences  assuming  the  coefficient  k^  depends  only  X2.  In  the  resulting  matrix 
operator  we  replace  the  Sommerfeld-like  boundary  condition  (4.8)  with  either  Dirichlet  or 
Neumann  boundary  condition.  The  resulting  matrix  is  then  used  as  a  pre-conditioner  to 
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accelerate  the  convergence  of  GMRES.  To  invert  this  matrix  effectively,  eigenvectors  of  the 
corresponding  operator  are  found  and  used  as  an  orthonormal  basis.  Such  an  orthonormal 
basis  is  possible  to  find,  because  both  Dirichlet  and  Neumann  boundary  conditions  generate 
self  adjoint  matrices,  which  is  the  key  reason  of  replacement  of  Sommerfeld  boundary  con¬ 
ditions  with  either  of  the  above  two  in  this  pre-conditioner.  This  leads  to  a  diagonal  matrix 
to  invert,  which  obviously  can  be  done  very  rapidly.  Iterations  are  used  "with  respect”  to 
both  boundary  conditions  (in  order  to  get  boundary  conditions  (4.8)  in  the  end),  and  the 
mine-like  targets.  In  order  to  generate  data  on  many  frequencies,  a  high  order  extrapolation 
technique  is  used  to  give  a  starting  value  for  a  higher  frequency  given  values  at  the  lower 
ones.  This,  in  combination  with  the  pre-conditioner,  leads  to  a  total  CPU  time  of  about  5 
minutes  for  150  frequencies  on  a  199  x  199  grid. 

The  idea  of  the  pre-conditioner  was  published  in  [3];  and  complete  details  were  published 
in  [1].  It  is  worthwhile  to  mention  here  that  this  work  received  a  warm  reception  from  the 
Editorial  Board  of  Journal  of  Computational  Physics  [1],  and  one  of  figures  of  [1]  was  selected 
for  the  cover  of  the  issue. 

In  all  numerical  examples  below  data  simulation  for  the  inverse  problems  were  provided 
by  solution  of  the  forward  problem  (4.5),  (4.8)  using  the  method  of  [1].  In  particular,  some 
results  of  solutions  of  this  problem  by  this  method  are  presented  on  Figs.  2-4. 

5.2  TASK  2. 

DEVELOPMENT  OF  A  SECOND  GENERATION  OF  THE  ESM  FOR  THE 
SOLUTION  OF  THE  INVERSE  PROBLEM 

We  refer  to  this  algorithm  as  the  ”p-method”  because  we  compare  it  below  with  an  analog 
of  Natterer’s  method,  which  we  call  ”H-method”.  Results  of  this  section  are  published  in 
[12]. 

As  it  is  always  the  case  in  the  field  of  inverse  problems,  it  is  useful  to  establish  an 
uniqueness  result  first.  This  result  is  formulated  in  Theorem  5.1. 

Theorem  5.1.  For  X2  >  0  let  =  uj^k^g,  where  the  complex  constant  kos  was  defined 
in  (4-11),  does  not  depend  on  u  and  x,  Ke{kos)  >  0  and  Im(A;os)  >  0.  Suppose  the  functions 
ip{xi,w)  and  ■0(a;i,a;)  in  (4-12a,h)  are  given  for  Xi  E  (—00,00)  and  u  E  (wmimt^max))  where 
0  <  Wmin  <  Wmax-  Then  there  exists  at  most  one  pair  of  functions  {v,  h)  satisfying  (4-6), 
(4-12),  (4-H)  in  R\  such  that  h{x)  is  a  bounded  function  with  compact  support  in  {x2  >  7} 
with  a  positive  constant  7,  and  v  E  iJ^’~*(i?+)  with  s  >  1/2. 

An  obvious  inconvenience  of  the  equation  (4.14)  is  that  this  is  one  equation  with  two 
unknown  functions  h{x)  and  v{x,u)  in  it.  The  key  idea  of  the  ESM  is  to  eliminate  the 
unknown  perturbation  term  h{x)  from  (4.14)  using  differentiation  with  respect  to  a  ’’free” 
parameter  u>.  This  leads  to  an  integro-differential  equation,  in  which  integrals  are  taken  with 
respect  to  u  (see  below),  supplied  by  the  boundary  condition  resulting  from  (4.8),  (4.12). 
Naturally,  the  next  question  is:  How  to  solve  the  resulting  boundary  value  problem  (BVP)? 

In  the  first  generation  of  this  ESM,  which  was  aimed  on  medical  applications  and  was 
applied  to  the  time  dependent  diffusion  equations,  integrals  were  eliminated  via  truncated 
generalized  Fourier  series  with  respect  to  time  [10-14].  As  a  result,  a  BVP  for  an  elliptic 
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system  was  obtained  (the  latter  led  to  the  name  of  this  method).  A  shortcoming  of  this  spe¬ 
cific  implementation  was  that  while  locations  of  abnormalities  were  imaged  very  accurately, 
values  of  unknown  coefficients  within  them  were  estimated  with  a  poor  accuracy.  This  is 
clearly  unacceptable  for  the  applications  to  imaging  of  land  mines,  in  which  an  important 
input  in  the  process  of  differentiation  of  land  mines  from  clutter  would  come  from  the  values 
of  electrical  parameters  within  those  targets.  This  motivated  the  development  of  the  sec¬ 
ond  generation  of  the  ESM,  in  which  the  BVP  for  the  integro-differential  equation  is  solved 
directly,  i.e.,  without  elimination  of  integrals.  The  resulting  algorithm  provides  accurate 
images  of  both  locations  of  targets  and  values  of  electrical  parameters  within  them. 

5.2.1.  mTEGRO-DIFFERENTIAL  EQUATION 


Note  that  because  of  (4.4),  the  function  uq  in  (4.14)  has  the  form 

Uo  =  T{ko)  exp{iukosX2)-  (5.1) 

The  function  h{x)  in  (4.14)  can  be  isolated  by  dividing  both  sides  by  uj^kl^uo.  Let 

uPk^^UQ 

Then  because  of  (5.2),  (4.14)  becomes 

V^H  +  2iujkosH^2  =  ^  e  (5.3) 


with  the  corresponding  boundary  conditions 

^\x2=o  ~  ‘P{Xl,(x)),  -^12112=0  ~  ^  ^  (^min)  ^<^inax))  (5.4a) 

±  ii^kosH\xi=±Li  =  0  (5.4b) 

^\x2=L2 

where  functions  (p  and  are  obtained  from  the  functions  p  and  ip  of  (4.12)  in  an  obvious 
way. 

To  eliminate  the  function  h{x)  from  (5.3),  differentiate  this  equation  with  respect  to  u. 
Let 


We  assume  that  for  every  x  € 
and 


II 

lim  H{x,  u)  =  0 

LJ—^OO 

(5.6) 

P^Px2  ^  Oo), 

(5.6) 
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as  functions  of  uj.  In  the  case  when  the  2-D  boundary  value  problem  (4.5),  (4.6)  is  consid¬ 
ered  over  the  whole  space  including  the  case  of  the  linearized  (with  respect  to  h[x)  ) 
equation  (4.14),  conditions  (5.5)  and  (5.6)  easily  follow  from  the  above  mentioned  result  of 
[3].  However,  at  this  point,  we  have  not  been  able  to  rigorously  establish  these  conditions 
for  the  case  of  the  finite  domain  fi.  Still,  we  have  observed  them  in  the  computations  of 
solutions  of  the  forward  problem.  We  assume  below  that  conditions  (5.5)  and  (5.6)  hold. 
Hence, 

OO 

H{x,lo)  =  —  J  p{x,T)dT. 

UJ 

However,  we  will  use  approximation 

OO  £^max 

p{x,T)dT  w  /  p{x,T)dT 

LJ  UJ 

and  thus 

Wmax 

H{x,u))  =  —  J  p{x,T)dT.  (5.7) 

UJ 

In  (5.7)  is  replaced  with  “=”  for  the  sake  of  notational  convenience. 

Let  giixi,uj)  =  ■^(p{xi,u)  and  g2{xi,(jj)  =  £;'4>{xi,u).  Then  because  of  (5.3),  (5.4)  and 
(5.7),  we  obtain 


tA^max 

/* 

V'^p+2i(jjkosPx2  =  2iA:os  J  Px2{^^'^)dr,  x  €  fi, 

(5.8a) 

UJ 

pU2=o  =  gi{xi,oj),  Px2l^=o  =  92{Xi,U}), 

(5.8b) 

UJxnBX 

/» 

Px^  T  ^^^05^1x1  =±1/2  J  p{,^j'^)dT\xx=:kLi'> 

(5.8c) 

UJ 

Px2\x2=^L2  ~ 

(5.8d) 

Thus,  we  have  obtained  the  boundary  value  problem  (5.8)  for  the  integro-differential  equation 
(5.8a)  with  Volterra-like  integrals  being  present  in  both  the  equation  itself  and  the  boundary 
conditions.  Once  the  function  p  is  found  from  (5.8),  the  function  h{x)  can  be  easily  recovered 
by  backwards  calculations:  First,  the  function  H{x,uj)  is  available  through  (5.7);  next,  the 
function  h{x)  can  be  calculated  from  (5.3)  evaluated  at  w  =  tUmin-  Therefore,  the  principal 
computational  question  becomes:  How  to  solve  the  boundary  value  problem  (5.8)7  This 
question  is  addressed  below. 

The  problem  (5.8)  is  overdetermined,  because  both  Dirichlet  and  Neumann  boundary 
conditions  are  given  at  0:2  =  0  in  (5.6).  There  is  no  guarantee  that  the  solution  of  the 
overdetermined  problem  (5.8)  exists.  Therefore,  the  idea  is  to  find  such  a  function  p,  which 
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satisfies  conditions  (5.8)  in  an  optimal  sense.  That  is,  we  want  to  find  a  minimal  norm 
solution.  In  terms  of  matrix  equations  resulting  from  the  finite  difference  approximation  of 
(5.8)  for  each  value  of  oj,  this  means  that  we  will  find  a  normal  solution  of  the  matrix  system 
for  each  oj,  which  is  a  parameter  of  this  system. 

5.2.2.  AN  IDEALIZED  CASE  OF  COMPLETE  DATA  COLLECTION 


To  explain  our  idea  better,  we  assume  in  this  subsection  that  both  Dirichlet  and  Neumann 
boundary  conditions  are  given  over  the  entire  boundary  d^l.  Then  (5.8)  becomes 


A‘^{p)  :=  +  2iu}kosPx2  =  2*fcos 


Wmax 

j  Px2{x,T)dT, 

UJ 


X  eQ, 


(5.9a) 


plan  =  gi{x,aj),  ~  92{x,uj), 


(5.9b) 


where  is  the  differential  operator  in  the  left  hand  side  of  (5.8a),  and  the  functions  g\  and 
P2  are  given  for  {x,u})  e  dQ  x  (wmiint^max),  with  pi|i2=o  =  9i,P2U2=o  =  92-  First,  consider 
the  more  general  problem 

A^{P)  =  S{x,uj),  (5.10a) 


•Flan  9l) 


9Pi  _  ~ 

^|an  —  92, 


(5.10b) 


where  the  function  S  E  1-2(0  x  (wmima^max))- 

Important  Remark.  We  note  that  if  the  right  hand  sides  of  equations  (5.3),  (5.4) 
are  zeros,  then  H  =  0.  Similar  statement  is  true  for  the  system  (5.8):  if  the  right  hand 
sides  of  equations  (5.8)  equal  zero,  then  p  =  0.  This  follows  from  the  well  known  theorem 
about  uniqueness  of  solution  of  the  so-called  ’’Cauchy  problem”  for  the  elliptic  equation, 
cf.  [15,  16].  This  observation  is  important  for  our  computations,  because  it  justifies  linear 
independence  of  column  of  matrices  A^  and  B{/3j)  below,  provided  that  the  grid  size  is 
sufficiently  small.  This,  in  turn,  implies  that  matrices  A2*A2  and  B*{Pj)B{Pj)  are  Hermitian 
Positive  Definite  ones,  which  enabled  us  to  use  Conjugate  Gradient  Method  for  their  solution. 

We  seek  solution  of  the  problem  (5.10)  in  the  minimal  norm  sense  as 


\\A‘^{P)  -  5'(x,a;)|(^^(jj)  ^  min,  for  every  a;  G  (  ^min,  a^max), 


F|afj  9l, 


du 


Ian  =  ^2) 


(5.11a) 

(5.11b) 


P{x,Lj)  E  for  each  fixed  u  E  [caminitamax]  and  sup  ||P||//2(n)  <  oo  (5.11c) 

[u-niin  j^max] 

This  minimization  problem  is  equivalent  to  the  solution  of  the  4th  order  elliptic  PDE; 

{A'^*A‘^){P)  =  A'^^S),  (5.12) 
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plus  the  conditions  (5.11b).  Next,  we  replace  in  (5.12)  the  function  S  with  the  right  hand 
side  of  (5.9a).  Thus,  we  finally  come  up  with  the  following  boundary  value  problem  for  an 
elliptic  integro-differential  equation  of  the  4th  order 


f  C*^max  ^ 

(P)  =  A'^*  2ikos  J  Px2{x,T)dT  ,  xgQ, 


dP, 


Plan  =  a),  a), 

plus  the  condition  (5.11c) 

The  problem  (5.13)  can  be  solved  iteratively  as 


(Pn)  =  A‘^* 


2ikos  j  i^Pn- 


u> 


L {^,T)dT 


(5.13a) 

(5.13b) 

(5.13c) 


(5.14a) 


Po  :=  0  (5.14b) 

9Pn 

Pn\dsi  =  9i{x,u}),  (5.14c) 

Therefore,  on  each  iterative  step  n  one  should  solve  the  boundary  value  problem  (5.14)  for 
the  4th  order  elliptic  operator  =  A‘^*A^  for  each  value  of  the  parameter  uj  G 

To  formulate  convergence  result  for  this  iterative  process,  we  first  introduce  the  Banach 
spaces  C\  =  C  ^Co'''^(fl),  ^^(fi)^  and  =  C  ^C'^(fl),  Co'^^(0)j  of  bounded  linear  opera¬ 
tors.  Let  =  A^*A^.  Obviously,  B^^  ^  C\.  The  existence  of  the  operator  B~^  €  £2  for 
every  u)  G  [wmimt^^max]  follows  from  Theorem  4.1  of  [11].  Furthermore,  the  perturbation 
theory  for  linear  operators  implies  that  sup  ||Pj^||  <  with  a  positive  constant  M, 

cf.  Theorem  2.23  in  Chapter  4  of  [15].  We  also  assume  that  there  exists  a  function  F{x,u>) 
such  that  _ 

1.  F{x,u)  G  for  every  uj  G  [cUminj^max])  and  the  norms  ||F(x,a;)||4+^  are  uni¬ 

formly  bounded  for  all  uj  G  [wmimn^max]-  Here  /?  =  const.  G  (0, 1). 

2.  PUn  =  h{x,uj)  and  =  g2{x,uj). 

A  method  of  construction  of  such  a  function  F  was  described  in  [13].  The  following 
Theorem  can  be  proven  using  the  above  result  concerning  the  operator  B~^  as  well  as  the 
conventional  approach  to  Volterra-like  integral  equations: 

Theorem  5.2.  Let  n  C  be  a  convex  bounded  domain  with  dQ,  G  C°°.  Also,  assume 
the  existence  of  the  function  F{x,  uj)  with  the  above  properties.  Further,  let  the  operators 
A'^  have  the  form  (5.9a)  and  B^  =  A^*AP.  Then  for  every  uj  G  [a;o,t«^i]  there  exists  an 
inverse  operator  B~^  G  and  sup  ||Pj^||M,  where  M  is  a  positive  constant;  the 
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boundary  value  problem(5.13)  has  an  unique  solution  P*  such  that  G  for  every 

u)  G  [c^oj^i];  the  iterative  procedure  (5.14)  converges  to  this  solution;  and 

sup  ||P|c  ””  PnlU+Z?  ^  (^max  ““  ^min)  j 

W€[Winin)t^max] 

where  the  positive  constant  K  depends  only  on  the  domain  Vt,  the  constant  kos  in  (5.9a)  and 
the  numbers 


5.2.3.  INCOMPLETE  DATA  COLLECTION 

Briefly,  the  idea  here  is  as  follows.  Let  Wmin  —  /3o  <  Pi  <■■■<  Pn  =  t^max  be  a 
discretization  of  the  frequency  band  [u^miiDt^max]  with  a  uniform  step  Au.  Let  :=  A^  be 
matrix  representing  the  finite  difference  analog  of  the  operator  +  2iu)kos^  in  (5.8a),  as 
well  as,  the  left  hand  sides  of  boundary  conditions  (5.8  b-d).  Let  p^  be  the  corresponding 
discrete  approximation  of  the  function  P{x,u})  at  a;  =  pn,  and  p  =  (p°,  •  •  -p^)-  Further,  let 
S'"(p)  be  the  matrix  representing  the  finite  difference  approximation  of  the  right  hand  sides 
of  (5.8  a-d),  in  which  integrals  are  taken  from  u  =  Pn  to  Wmax  in  a  discrete  form.  Then  the 
discrete  analog  of  the  system  (5.8)  is  the  following  matrix  system 

AW)  =  S2{p)  (5.15) 

We  assume  that  values  p^  for  A:  =  n  +  1,  •  •  •  ,  A  are  given  and  p^  =  0.  The  system  (5.15) 
is  overdetermined,  because  of  two  boundary  conditions  (5.8b),  rather  than  one.  This  means 
that  the  number  of  rows  in  the  matrix  A'l  exceeds  the  number  of  columns.  Thus,  we  seek  a 
normal  solution  of  this  system  as 


{ATA-,){r)=^Arsm 


This  was  a  brief  description  of  our  idea  for  the  case  of  incomplete  data  collection.  Details 
are  given  below  in  this  section. 

Consider  the  second  order  central  finite-difference  approximation  of  the  boundary-value 
problem  (5.8)  with  an  uniform  mesh  cell  size  of  x  where  h^^  =  hx^  =  and 
Mxi  and  Mx2  are  the  number  of  grid  points  in  the  Xi  and  X2  directions,  respectively.  The 
gridded  region  uses  Xi  values  over  the  interval  [— Li  +  hxi/2,  L\  —  hij/2]  and  X2  values  over 
[/ii,2/2,I/2  —  hx2/2].  The  boundary  conditions  (4.8b)-(4.8d)  are  imposed  by  use  of  a  second 
order  approximation  formula  centered  on  each  boundary,  using  fictitious  values  outside  of 
fi,  which  are  eliminated  in  setting  up  the  matrix  system.  In  this  case  the  discretized  system 
(5.0)  can  be  written  in  the  following  form 

Tf  =  /  A.,p(r)<iT  + 

CJ  Uf 
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where  A^,  is  the  central  second  order  approximation  of  the  first  derivative  in  the  X2  direction; 


where  ( 


where  /  ’ 

\  df.  =  0,  otherwise 


^21 J  •  *  *  ) 

dij ) . . . 

■  ) 

=  1,- 

1,- 

• . ,  Mxi , 

11 

=  1 ,  .  .  .  ,  Mx2 

'rfu, 

^21 5 • • *  ? 

4,... 

■  > 

Mx^  J 

)  3 

=  2,- 

*  •  J  ^X2  . 

,  g\f^  ^  and  g'^  =  (^gl, . . . ,  glf^  j  are  approximations  of  the  boundary  condi¬ 
tions  (5.8b)  and G  =  (Gn, . .  • , i, G12,  •  •  • , Gmx^2, 0 . . . , 0)^,  where Gn  =  gj  —  -^g\,Gi2  = 
{iashx^  -  1)2^^ 

In  (5.16),  is  the  rectangular  matrix,  which  consists  of  the  finite  difference  approxima¬ 
tions  of  the  left  hand  sides  of  equalities  (5.8a-d),  multiplied  by  The  number  of  rows  in 
the  matrix  exceeds  the  number  of  columns,  because  of  the  overdetermination  in  (5.8b). 
To  find  the  standard  least  squares  solution  of  the  overdetermined  system  (5.16),  we  use  the 
method  of  normal  equations  as  in  equation  (5.18)  below.  In  practice  the  final  normal  system 
is  similar  to  the  system,  which  would  result  from  setting  up  the  4th  order  finite  difference 
operator  as  in  (5.13a).  However  the  advantage  of  our  treatment  of  the  incomplete  data  case 
is  that  no  explicit  treatment  of  an  unknown  second  boundary  condition  along  parts  of  the 
boundary,  other  than  {x2  =  0},  is  required.  In  particular  no  spatial  weight  function  needs 
to  be  introduced,  as  it  was  the  case  in  earlier  works  [13,15,16] 

In  principle,  the  system  (5.16)  can  be  solved  by  iterations,  similar  to  (5.14).  But  because 
it  is  Volterra-like  integro-differential  equation,  a  more  natural  way  is  to  use  the  marching 
method  [4],  which  produces  the  solution  in  one  step.  First  we  approximate  the  integral 
terms  in  the  right  hand  side  of  (5.16)  by  a  simple  trapezoidal  rule.  For  this  purpose  we  use 
a  regular  mesh  in  the  ^-direction: 


Pu  =  (^min  +  where  n  =  0, . . .  ,N  and  A cu  = 


Wn 


- 


N 


Then  (5.16)  can  be  rewritten  in  the  following  form 

iV-l 


_  2ikshl^D^  A; 


W+P 


0^+1 


X2 


•  Ao;  + 


Aik. 


2  (2  i'Pnkshxi)'^  hxy 

%~n 

Finally,  we  obtain  the  system 

/ly  =  s;(p"+‘. 


P  +  P . A  a;  -f  G". 


(5.17) 
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where 


A2  =  A^"-iKhlAuD^A,, 


and 


2,ikg  ^X2 

(2  —  Pn^s^xi)"^  hxi 


AcoD^ 


N-l 

S2  =  2ishl^AuD^  Y,  + 

i—nA-l 


4iks 

(2  i^n.kshxi')^  hxi 


N-l 

AuD^  p  +  G^, 

i=n+l 


taking  into  account  the  assumption  that  =  0.  Note  that  the  matrix  A^  depends  only  on 
the  frequency  Pn,  but  not  on  the  right  hand  side,  or  the  solution  p  for  any  value  of  j. 

Similarly  to  (5.16),  (5.17)  is  an  overdetermined  system.  To  solve  it  for  each  starting 
from  n  =  N  —  1  to  n  =  1,  we  use  the  normal  equations  method.  Thus  [  ],  we  replace  (5.17) 
with  the  normal  equations  form 


ATA'S(r)  =  ATS2(P).  n  =  Ar-l,...,l.  (5.18) 

Thus,  on  each  step  n  of  the  marching  method  we  must  solve  the  system  (5.18). 

Iterations  with  Respect  to  H{x,Ujaax) 

The  above  described  procedure  relies  on  the  formula  (5.7),  which  is  an  approximation  for 


H{x,U})  =  -  J  p{x,T)dT  +  H{x,U}inax)-  (5.19) 

U) 

So,  the  use  of  (5.7)  implies  that  the  term  H{x,uJmax)  in  (5.19)  is  neglected.  In  our 
numerical  experiments,  we  first  tested  the  case  where  this  term  was  included.  Namely,  we 
used  an  accurate  value  for  H{x,ujjasx),  which  was  computed  by  the  solution  of  the  above 
forward  problem  of  data  simulation.  In  such  a  case  the  above  described  inverse  algorithm 
led  to  images,  which  were  almost  identical  to  the  correct  ones.  Therefore,  we  decided  to 
mitigate  the  influence  of  the  cut-off  of  the  function  H{x,Ujnax)  by  the  following  iterative 
procedure. 

Step  1.  We  use  formula  (5.7),  i.e.,  we  assign  H{x,uimax)  ■=  0  and  compute  the  function 
Pi{x,u})  :=  P{x,u!)  for  u  E  by  the  procedure  described  above.  Given  Pi{x,a), 

we  compute  Hi{x,a)  as 


Hi{x,u)) 


ttJmax 

J"  ■fl(^)^)d'r,  U  E  )  ^max] 

u> 


Next,  given  the  function  Hi{x,  u),  we  compute  the  perturbation  term  hi{x)  using  Hi{x,oj) 
for  the  lowest  value  of  uj Wmin,  that  is  by  (5.3) 


hi{x)  =  —  +  2iujksH\x2 


(5.20) 
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Finally,  given  the  function  hi{x)  from  (5.20),  we  solve  the  forward  problem  (4.5),  (4.7)  for 
the  function  v{x,u)^^.  Having  the  resulting  function  V\  :=  u,  we  compute  the  function 
Hi{x,u>iogx)  by  the  formula  (5.2). 

Step  n  >  1.  In  the  above  procedure  for  the  solution  of  the  inverse  problem,  we  modify 
(5.7)  as 

Hk{x,U))  =  -  j  Pk{x,T)dr  +  Hk-i{x,Umax). 

LJ 

Then  we  proceed  as  in  step  1. 

We  found  in  our  numerical  experiments  that  this  iterative  procedure  often  enables  one 
to  improve  the  quality  of  the  initial  image  of  h{x)  :=  hi{x). 

Another  important  problem  to  address  here  is  a  choice  of  a  fast  numerical  method  for 
the  solution  of  the  system  (5.18).  This  system  must  be  solved  for  many  frequencies  involved. 
In  addition  because  of  high  frequencies  involved,  we  need  to  use  a  fine  spatial  grid,  which  is 
similar  to  the  above  discussed  case  of  the  forward  problem.  Thus,  if  conventional  Gaussian 
elimination  techniques  were  used,  this  would  be  a  time  consuming  procedure. 

The  first  key  issue  to  resolve  here  is  the  band  limited  structure  of  the  matrix  = 
A2*A2,  since  it  is  generated  by  a  differential  operator  (rather  than  an  integral  operator 
in  a  conventional  setting  of  inverse  problems).  The  second  key  issue  is  the  fact  that  the 
columns  of  A^  are  linearly  independent,  because  of  the  Important  Remark  above.  Hence  the 
matrix  C2  is  Hermitian  Positive  Definite.  Thus,  we  have  chosen  a  preconditioned  Conjugate 
Gradient  Method  using  the  method  of  nested  dissection  [8],  but  only  for  a  small  number  of 
frequencies.  In  this  approach  we  have  developed  an  automatic  algorithm  for  the  near  optional 
choice  of  frequencies  ranges,  over  which  we  use  the  same  preconditioner.  This  method  has 
enabled  us  to  dramatically  reduce  the  overall  computational  time.  For  a  discretization  of 
200x70,  using  132  frequencies  in  the  marching  method,  this  approach  yields  a  solution  of 
less  than  3  minutes  for  each  iteration  of  H{x,uJraax)  (usually  two  such  iterations  are  used) 
using  one  processor  on  a  SGI  Origin  200.  If  using  fast  dual  procesing  techniques,  this  timing 
can  be  decreased  by  a  factor  of  five.  But  even  the  above  timing  is  quite  acceptable  for  the 
above  diagnostic  procedure  of  land  mines.  More  details  about  this  method  can  be  found  in 
[2  ].  Algebraic  system  of  section  5.5  is  also  solved  by  this  method. 

5.3.  TASK  3. 

DEVELOPMENT  OF  THE  MODIFIED  VERSION  OF  THE  METHOD  OF 

TASK  2  AND  COMPARISON 
5.3.1.  Description  of  the  ff-method 

The  H-method  is  a  second  version  of  the  above  algorithm.  Unlike  the  first  version,  we 
do  not  use  the  differentiation  of  the  data  with  respect  to  the  frequency  u)  here.  This  seems 
to  be  more  suitable  in  the  view  of  possible  work  with  the  experimental  data.  In  addition, 
unlike  the  p-method,  the  H-method  does  not  use  an  assumption  that  the  loss  tangent  is 
independent  on  frequency.  Another  interesting  feature  of  the  H-  method  is  that  it  is  a 
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version  of  Natterer’s  method,  which  is  the  best  competing  technique  in  the  field  [5,6,19]. 
This  enabled  us  to  make  a  comparison  with  the  best  method  in  the  field. 

Let  h  be  the  grid  function  approximating  h{x)  and  B{u))  be  the  finite  difference  based 
matrix  representing  the  left  hand  sides  of  equations  (5.3),  (5.4)  as  a  function  of  oj.  I^s  iin- 
portant  to  note  that  B{u})  depends  only  on  oj,  kos  and  the  spatial  grid.  Also,  let  S^{u,  h,(p,  ip) 
be  the  finite  difference  operator  for  the  right  hand  sides  of  (5.3),  (5.4).  depends  not  only 
on  a>,  Uqs  and  the  spatial  grid,  but  also  on  a  grid  function  h  and,  of  course  on  the  overde¬ 
termined  boundary  data  (p  and  ip  evaluated  at  the  interface  grid  points.  Then  on  step  j  we 
solve  the  overdetermined  linear  system 

B{Pj)W  =  (5.21) 

where  the  vector  represents  the  finite  difference  approximation  to  the  grid  values  of 
H{x,l3j),  given  h.  Note  that  the  system  (5.3),  (5.4))  gives  one  equation  for  each  spatial  grid 
point,  plus  one  additional  equation  for  each  grid  point  along  the  air-soil  interface  X2  =  0. 

To  solve  the  system  (5.21),  we  use  method  of  normal  solutions,  thus  coming  up  with  the 
system  _ 

B*{0j)B{Pj)H^  =  B*{Pj)S^{(3j-,h,p,iP),  (5.22) 

The  system  (5.22)  is  solved  by  the  same  method  of  the  above  system  (5.18).  Note  that 
Hi  ^  i.e.,  systems  (5.21)  and  (5.22)  are  not  equivalent.  Next,  we  use  (5.1)  to  update  h 

as  _ 

h  :=  U=/3,  (5.23) 

where  the  differential  operators  are  understood  in  terms  of  finite  difference. 

The  complete  algorithm  of  the  i?-method  can  now  be  stated  as  follows:  Initialize  h  to 
be  the  zero  vector.  Then  complete  a  series  of  sweeps  (beginning  each  new  sweep  with  the 
latest  value  of  h  from  the  previous  sweep)  as  follows  until  the  stopping  criteria  is  satisfied. 

For  j  =  1 ..  .n  :  _ 

I.  Solve  the  system  (5.22)  using  the  latest  value  of  h. 

II.  Use  (5.23)  to  compute  the  updated  value  for  h. 

III.  Repeat  I-III  A:  >  1  times. 

5.3.2.  Brief  description  of  Natterer’s  method 

In  1995  F.  Natterer  and  F.  Wuebbeling  [19]  proposed  a  new  elegant  method  for  solutions 
of  inverse  problems  of  ultrasound  imaging.  This  method  is  also  based  on  resulting  differ¬ 
ential,  rather  than  conventional  integral  operators.  Later  this  idea  was  extended  on  other 
applications,  by  O.  Dorn,  et.al.  [5,6] . 

Briefly,  the  idea  of  the  method  [19]  is  as  follows.  First,  the  solution  of  an  overdetermined 
BVP  (similar  to  one  of  (5.3),  (5.4))  is  calculated  for  one  value  of  a  ’’free”  parameter  under 
an  assumption  that  the  perturbation  term  (like  h{x))  is  known.  This  free  parameter  is 
the  source  position  in  the  above  referenced  works,  and  is  frequency  cj  in  our  case.  That 
overdetermined  problem  is  solved  by  the  method  of  normal  solutions,  which  is  similar  to 
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the  above.  Given  the  normal  solution  of  this  BVP,  the  value  of  the  perturbation  term  is 
updated  as  the  right  hand  side  of  an  equation,  which  is  similar  to  (5.3).  Next,  this  updated 
perturbation  term  is  used  as  an  input  for  the  same  procedure  being  repeated  for  the  next 
value  of  the  free  parameter.  Thus,  iterations  are  performed  with  respect  to  the  value  of  the 
free  parameter.  It  was  shown  numerically  in  [5,6,19]  that  this  process  often  converges,  if 
the  unknown  perturbation  term  h{x)  is  sufficiently  small. 

This  process  was  naturally  called  in  [19]  ’’propagation  -  backpropagation  method”.  In¬ 
deed  first,  one  ’’propagates”  a  ’’pseudo”  field  into  the  medium  by  solving  the  overdetermined 
BVP  for  an  assumed  value  of  h{x).  Next,  one  ’’propagates  if  backwards”  by  updating  h{x). 

The  key  advantage  of  the  algorithm  [19]  over  conventional  techniques  is  that  it  solves 
differential,  rather  than  integral  equations  on  each  step,  which  is  similar  to  the  idea  of  the 
ESM.  The  major  disadvantage  of  the  specific  implementation  of  this  method  in  [5,6,19], 
however,  is  that  instead  of  using  the  entire  resulting  matrix  for  inversion,  only  diagonal 
elements  of  this  matrix  are  counted  in  these  publications.  This  would  mean,  that  in  the 
above  case  of  if-method,  only  diagonal  elements  of  the  matrix  A2*A2  in  (5.18)  would  be 
counted,  which  would  dramatically  decrease  the  quality  of  the  images.  Thus,  in  order  to  have 
a  fare  comparison,  we  decided  to  use  an  advanced  version  of  the  method  [19]  by  counting 
the  entire  resulting  matrix,  rather  than  its  part.  Although  we  are  not  making  a  direct 
comparison,  but  the  improvement  seems  clear. 

5.4.  NUMERICAL  EXPERIMENTS 

The  main  goal  of  the  numerical  experiments  presented  in  this  section  is  to  demonstrate 
and  compare  the  properties  and  performances  of  both  p  and  H  methods  for  realistic  ranges  of 
parameters  and  frequencies.  The  values  of  the  coefficients  in  the  Helmholtz  equation,  which 
correspond  to  the  electromagnetic  properties  of  air,  soil  and  different  targets  were  presented 
in  Table  I.  In  the  numerical  experiments  the  background  medium  is  a  wet  soil  with  a  5% 
moisture  content.  The  targets  are  assumed  to  be  filled  with  TNT. 

The  physical  domain  for  the  inverse  problem  is  selected  to  be  [—Lxi ,  Lij]  ^  [0)  -^*2]  where 
I/xj  =  60  cm  and  Lx^  =  40  cm,  with  a  201  x  71  grid.  In  both  approaches  the  system  was 
non-dimensionalized  in  space  and  frequency.  However  for  simplicity  the  results  are  presented 
here  in  the  original  coordinate  systems.  The  spatial  grids  selected  are  uniform,  but  do  not 
necessarily  have  the  same  spacings  in  xi  and  x^- 

5.4.1.  The  most  practically  important  case  of  a  simple  target  filled  with  TNT. 

First  we  consider  a  case  of  a  simple  mine-like  inclusion,  which  seems  to  be  the  most 
important  one  for  applications  to  diagnostics  of  mine-like  targets.  This  is  a  circular  target 
filled  with  TNT,  with  the  center  at  re  =  (rri,rr2)  =  (10,5)  cm  and  a  diameter  of  5  cm.  Thus 
h{x)  =  {klj^T  -  K>etsoii)lK,etsoii  =  -0-319  -  0.152i  The  real  part  of  the  corresponding 
function  h{x)  is  displayed  in  Figure  1. 

For  this  application  simulations  of  the  solution  of  the  forward  problem  are  made  using 
the  forward  solution  for  the  discrete  frequencies  uij  G  27r  •  [0.5,3]  GHz  with  the  step  size 
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Target  1:TNTin  wet  soil 


Boundary  cofxfitions  for  H  at  die  air-soil  interface 


Figure  1:  Re[/i(a:)]  for  a  circular  Figure  2:  Re[if(x,ct;)]  as  a  function 

shaped  signle  target  filled  with  TNT  of  the  frequency  /  at  x  =  (10, 0). 


of  27r  •  0.01  GHz.  The  physical  domain  is  taken  to  be  a  square  with  300  cm  sides  centered 
about  the  air-soil  interface,  and  a  400  x  400  point  computational  grid  is  used  to  achieve 
accuracy  at  the  higher  frequency  values.  Two  of  the  studies  to  be  made  here  consider  the 
effect  of  various  choices  of  w„,in,‘*^max  and  the  frequency  spacing  Aw  on  the  quality  of  the 
inverse  problem  solution. 


For  Figures  2,  3  and  4  the  solid  line  represents  the  real  parts  of  the  original  values  obtained 
through  the  solution  of  the  forward  problem,  the  stars,  which  are  the  data  originally  received 
by  the  algorithm,  represent  those  values  with  a  =  0.1  multiplicative  Gaussian  noise  added, 
and  the  dot-dash  line  the  result  after  a  cubic  spline  smoothing  process  [1].  Figures_2  and  3 
display  the  real  parts  of  H{x,u)  =  ^(10,  w)  and  the  normal  derivative  Hx2{x,u})  =  ■0(10,  w) 
as  a  function  of  f  =  —,  just  above  the  target  at  x  =  (10,0)  on  the  air-soil  interface. 

Despite  the  scatter  in  the  noisy  data,  the  cubic  spline  smoothing  process  appears  to  be 
doing  a  good  job  of  returning  smoothed  values  close  to  the  original.  The  largest  differences 
are  around  the  peaks  of  the  curves  and  to  the  far  left  of  Figure  3.  Figure  4  displays  the 
real  part  of  H{x,uio)  —  0(x,wo)  along  the  air-soil  interface  {x2  =  0}  for  wq  =  1.0  •  27r 
GHz.  It  should  be  clarified,  that  smoothing  was  done  for  each  spatial  point  x  =  xu,  where 
Xu  =  —Lxi  +  Aii(i  —  1/2), z  =  1.  •  •  •  ,Mxi,  with  respect  to  the  frequency  w  as  described 
above.  However,  no  smoothing  was  performed  with  respect  to  the  spatial  variable  Xj. 


In  this  first  test  the  p  and  H  imaging  algorithms  were  applied  to  a  single  target  of 
circular  shape  filled  with  TNT  as  specified  above.  The  results  used  data  over  frequencies 
from  0.5  GHz  to  3  GHz,  using  an  increment  of  A/  =  .01  GHz.  The  stopping  criteria 
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Figure  3:  Re[i/a:2(^)‘^)]  as  a  function 
of  the  frequency  /  at  x  =  (10,0). 


Figure  4:  Re[/f(xi,£u)]  for  uq  =  1.0  • 
2ttGHz  as  a  function  of  x\  along 
{X2  =  0}. 


Figure  5:  Re[himaged{x)]  for  the  solu¬ 
tion  of  the  inverse  problem  by  the  p 
method. 


Figure  6:  Re[/iimoped(a;)]  for  the  solu¬ 
tion  of  the  inverse  problem  by  the  H 
method. 


Horizontai  cross  section  of  h  (p  mettwd] 


Vertical  cross  section  of  h  (H  method} 


Figure  7:  Horizontal  cross-section  of 
^'^[Kmaged{x)]  after  two  sweeps  using 
the  p  method  for  different  values  of 


Figure  8:  Horizontal  cross-section  of 
Re[/iimafled(3^)]  after  two  sweeps  using 
the  H  method  for  different  values  of 

^max- 


here  and  elsewhere  requires  running  one  more  sweep  than  was  used,  stopping  when  the  last 
result  either  decreased,  showed  significant  oscillations  or  changed  very  little.  Figures  5  and  6 
display  the  real  part  of  the  imaged  function  h  :=  himagedi^a)  obtained  after  two  sweeps  of  the 
p-method  and  one  sweep  of  the  i7-method  respectively.  The  contour  plots  of  the  recovered 
function  Re(/i)  are  similar,  both  accurate  as  to  the  centered  location,  both  lacking  significant 
artifacts.  However,  the  maximal  value  of  the  H  solution  within  the  target  is  about  20%  off  the 
correct  one.  Whereas  the  value  for  the  p  solution  is  exact.  This  difference  speaks  favorably 
for  p  solution  in  the  light  of  the  goal  of  diagnostics  of  mine-like  targets.  The  imaginary  part 
of  h  here  and  elsewhere  is  in  general  less  satisfactory  because  |Im(fc)|  <<  |  Re(A;)|  by  Table 
1. 

The  effects  of  using  frequencies  from  fm\n  =  0.5  GHz  to  various  upper  values  of  /max  from 
1.0  GHz  to  3.0  GHz  was  also  considered  for  both  the  H  and  p  methods.  The  frequency 
spacing  A/  is  fixed  at  0.01  GHz.  To  clearly  demonstrate  the  results  quantitatively,  cross- 
sections  of  the  real  part  of  the  imaged  function  h  along  both  vertical  and  horizontal  lines 
are  displayed  for  eaoh  method.  In  each  case,  the  cross  section  is  along  the  line  where  the 
values  are  greatest.  Typically  for  the  horizontal  lines  this  is  close  to  X2  =  5  cm,  and  for  the 
vertical  lines  close  to  xi  =  10  cm.  In  Figures  7-10  the  solid  lines  represent  the  exact  values 
and  various  other  lines,  as  identified  by  the  legend,  the  cross  sections  of  various  computed 
solutions. 

As  one  can  see,  the  increase  of  the  range  of  the  frequencies  from  0.5  —  1  GHz  to  0.5  —  1.5 
GHz  to  0.5  —  2  GHz  gives  significant  improvements  in  the  heights  of  the  recovered  images 
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HorizonUP  cross  section  of  h  (p  methocl) 


y{cm) 

Figure  9:  Vertical  cross-section  of 
Re[/iiTOage(i(a;)]  after  two  sweeps  using 
the  p  method  for  different  values  of 

^nicix* 


Horizontal  cross  section  olh  (H  method) 


y(cm) 

Figure  10:  Vertical  cross-section  of 
Re[/ii^ased(3^)]  after  two  sweeps  using 
the  H  method  for  different  values  of 

^max- 


for  both  the  p-  and  if-methods  and  is  thus  clearly  desirable. 

The  effect  of  various  frequency  spacing  A/  =  Auj/2tt  was  also  examined.  It  was  con¬ 
cluded  that  it  is  unnecessary  to  use  a  fine  frequency  spacing  of  A/  =  0.01  GHz,  since 
A/  =  0.02  GHz  provided  the  same  results. 

5.4.2  Three  targets  filled  with  TNT 

In  the  remaining  numerical  examples  the  application  of  both  algorithms  to  the  case  of 
three  multiple  mine-like  targets  of  different  sizes  and  soil  depths  is  considered.  These  targets 
are  again  in  wet  soil  with  a  5%  moisture  content  and  filled  with  TNT,  using  parameter 
values  from  Table  I.  The  buried  objects  chosen  for  this  test  are  three  rectangular  mine-like 
targets.  Two  of  the  targets  are  5  x  4  cm  and  the  third  is  10  x  4  cm.  Three  mine-like  targets 
were  examined  to  see  if  the  H  and  p  inversion  algorithms  could  separate  multiple  scatterers 
and  reconstruct  well  the  deeper  object.  The  frequency  range  in  this  test  is  from  0.5  to  3.0 
GHz  for  the  p  method  and  from  0.5  to  2.0  GHz  for  the  H  method,  and  the  frequency  step 
is  A/  =  0.02  GHz  as  suggested  above.  The  two  smaller  targets  are  centered  5  cm  into 
the  ground  and  the  larger  rectangular  target  is  centered  10  cm  deep  into  the  ground.  The 
horizontal  centers  are  at  —10, 0  and  10  cm.  As  in  the  first  example,  the  detector  readings  are 
simulated  from  the  forward  problem  with  the  addition  oi  a  =  .10  multiplicative  Gaussian 
noise. 

The  real  part  of  the  corresponding  function  h{x)  is  displayed  in  Figure  11.  The  recon¬ 
structed  images  of  the  real  part  of  the  coefficient  are  shown  in  Figures  12-14.  From  these 
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Real  part  of  exact  image 


x(cm) 

Figure  11:  Re[/i(x)]  for  three  mine-like  targets  of  various  sizes. 

figures  it  can  be  seen  that  both  algorithms  perform  reasonably  well.  The  locations  and 
shapes  of  the  objects  are  both  fairly  accurate.  The  methods  provide  maximal  values  of  Re[h 
imaged  (2:)]  within  targets  which  are  about  12%  off  the  target  value  for  p-method  and  7%  off 
for  if-method. 

5.4.3  Conclusions  for  Comparison  of  p  and  H  Methods. 

The  iJ-method  generally  works  better  for  the  lower  upper  frequency  than  the  p-method. 
Also,  unlike  p-method,  i^-method  does  not  require  neither  the  differentiation  of  the  data 
with  respect  to  the  frequency,  nor  independence  of  the  loss  tangent  from  the  frequency, 
which  might  be  an  advantage  for  the  case  of  an  experimental  data.  Both  methods  provide 
accurate  locations  of  targets.  Another  key  parameter  for  the  goal  of  diagnostics  of  mine¬ 
like  targets  is  the  maximal  value  of  Re[h  imaged  (a:)]  within  targets.  In  the  case  of  a  single 
target,  which  is  the  most  important  one  for  a  practical  scenario,  the  p-method  provides  very 
accurate  maximal  values  of  Re[h  imaged  (a:)]  within  a  target,  as  opposed  to  the  iJ-method, 
for  which  those  values  are  20%  off  the  correct  ones.  The  latter  speaks  favorably  for  the 
p-method  in  the  light  of  the  above  diagnostics  goal. 

As  to  the  comparison  with  the  best  algorithm  in  the  field,  iif-method  is  an  advanced 
version  of  the  original  algorithm  of  F.  Natterer  [5,6,19].  In  those  works  only  diagonal 
elements  of  resulting  matrices  to  be  inverted  were  counted,  as  opposed  to  our  case,  in  which 
the  entire  matrix  is  counted.  Therefore,  if  we  would  literally  follow  the  idea  of  Natterer 
by  counting  only  diagonal  elements  of  those  matrices  in  the  above  H  method,  then  images 
would  be  much  worse  than  those  obtained  above. 

5.5.  CONCLUSIONS 
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Figure  12:  Final  reconstruction  of 
'RB[hijnaged{x)\  after  One  sweep  using 
the  p  method. 


Figure  13:  Final  reconstruction  of 
Re[^tma9ed(a^)]  after  one  sweep  using 
the  H  method. 


The  problem  of  differentiating  mines  from  clutter  using  hand-held  GPR  was  considered  as 
an  inverse  problem.  In  this  approach  locations  and  electrical  properties  of  mine-like  targets 
would  be  imaged  via  solution  of  an  inverse  problem.  The  input  data  for  such  a  problem 
would  be  the  back-reflected  electrical  signal  measured  at  many  frequencies  ranging  from  0.5 
to  3  GHz.  Locations  and  electrical  properties  of  targets,  in  turn  might  be  used  on  a  later 
stage  as  an  input  for  a  procedure  of  classification  of  targets. 

Three  major  tasks  were  achieved  in  this  project.  Each  of  these  required  about  one  year 
of  effort.  First,  a  new  rapid  algorithm  for  the  solution  of  the  forward  problem  on  high 
frequencies  was  developed  and  implemented  [1,3].  This  task,  though  axillary,  was  necessary 
to  perform,  because  the  conventional  classical  algorithm  of  Gaussian  elimination  is  too  slow 
for  this  range  of  frequencies. 

Next,  a  second  generation  of  the  Elliptic  Systems  Method  (ESM),  which  was  used  earlier 
in  Diffusion  Tomography,  was  developed  and  tested  for  the  inverse  problem  of  diagnostics  of 
mine-like  targets  using  GPR  [2]  (p-method).  The  key  new  element  of  this  second  generation 
method  is  that  the  resulting  integro-differential  equation  is  solved  directly,  rather  than  by 
eliminating  integrals,  as  it  previously  was.  In  addition,  a  fast  preconditioner  was  developed 
for  the  Conjugate  Gradient  Method,  which  significantly  speeded  up  solutions  of  inverse 
problems  of  this  project.  The  resulting  algorithm  provides  accurate  images  of  both  locations 
of  targets  and  values  of  the  Re[/i  imaged  (s;)]  within  them  in  about  6  minutes  of  CPU  time 
on  a  SGI  Origin  200  with  a  single  processor.  This  timing  is  well  acceptable  for  the  goal 
of  diagnostics  of  mine-like  targets,  since  the  process  of  diagnostics  is  much  slower  than  one 
of  screening.  Still,  this  timing  can  well  be  decreased  by  a  factor  of  five,  if  using  fast  dual 
precessing  techniques. 
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Horizontal  cross  section  of  h 


Figure  14:  Horizontal  reconstruction  of  Re[/ii„,aged(3j)]  using  the  p  and  H  methods. 

Finally,  a  modified  version  of  the  algorithm  [2]  was  developed  (ff-method),  which  is 
actually  an  advanced  variant  of  the  method  of  a  well  known  German  scientist  F.  Natterer 
[19].  It  was  shown  that  both  versions  provide  the  same  accurate  locations  of  targets,  while 
the  second  version  usually  performs  better  for  the  lower  values  of  the  upper  frequency. 
However,  in  the  practically  most  important  case  of  a  single  target  the  first  version  performs 
much  better  for  the  key  parameter  Re[h  imaged  (a;)]  within  the  target.  It  is  also  reasonable 
to  conclude  that  if  using  the  original  method  of  [19],  in  which  only  diagonal  elements  of 
matrices  are  counted  (as  opposed  to  its  more  advanced  version  used  here),  the  quality  of 
the  resulting  images  would  be  much  worse  than  ones  presented  above.  This  speaks  in  favor 
of  the  second  generation  of  the  ESM,  as  compared  with  the  best  competing  technique  of  F. 
Natterer. 


5.6  POSSIBLE  FUTURE  DIRECTION  OF  RESEARCH 

First,  it  would  be  very  interesting  to  verify  the  performance  to  both  p-  and  ^/-methods  on 
the  experimental  data.  The  methods  are  now  mature,  although  a  few  new  features  should  be 
implemented  before  this  would  be  possible.  Some  of  these  features  are:  the  measurement  line 
should  be  “raised”  from  the  air/ground  interface,  geometric  irregularities  of  the  air/ground 
interface  need  to  be  incorporated  in  the  model,  and  possibly  a  point  source,  as  opposed  to 
the  current  initializing  plane  wave,  should  be  incorporated. 

There  is  also  a  second  pressing  need  of  research,  connected  with  the  development  of 
globally  convergent,  as  opposed  to  locally  convergent,  inverse  algorithms  for  imaging  of  land 
targets.  This  would  be  a  far  advanced  second  generation  of  the  above  methods.  The  vast 
majority  of  inverse  techniques,  including  the  ones  presented  above,  is  heavily  relying  on  an 
assumption  that  the  properties  of  soil  are  known  with  a  good  accuracy.  This  leads  either 
to  the  linearization  of  the  inverse  problem,  or  to  a  perturbation  approach,  which  is  a  slight 
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modification  of  the  linearization.  In  practice,  however,  electrical  properties  of  soil  are  not 
always  known  with  a  good  accuracy.  Therefore  linearization  of  an  inverse  problem  is  not 
feasible  in  such  cases.  The  key  drawback  of  the  conventional  least  squares  objective  functions 
is  that  they  suffer  from  the  problem  of  local  minima,  as  soon  as  the  starting  vector  is  chosen 
far  from  the  solution  [18].  Therefore,  in  the  scenarios  when  properties  of  soil  are  not  known 
with  a  good  accuracy,  there  is  a  pressing  need  for  algorithms,  which  would  neither  rely  on 
the  perturbation  approach,  nor  suffer  from  the  problem  of  local  minima. 

A  new  idea  of  globally  convergent  algorithms  has  been  developed  by  M.  V.  Klibanov 
in  the  past  several  years,  cf.  [17].  It  was  only  recently  however,  when  the  computational 
feasibility  of  this  idea  was  verified  by  M.  V.  Klibanov  and  A.  Timonov  [18].  In  [18],  this 
approach  was  termed  convexification.  By  the  convexification  approach,  the  original  inverse 
problem  is  reduced  first  to  a  boundary  value  problem  for  a  non-linear  integro  differential 
equation,  in  which  the  unknown  coefficient  is  not  involved.  This  equation  is  similar  to  the 
above  equation  (5.8a),  except  of  the  non-linearity,  which  is  obviously  not  present  in  (5.8a). 
Next,  this  boundary  value  problem  is  solved  via  a  least  squares  cost  functional  J.  The 
key  new  element  of  J,  however,  is  the  presence  of  the  so-called  Carleman  weight  function 
(CWF).  The  role  of  the  CWF  is  that  it  guarantees  strict  convexity  of  J  on  a  compact  set  of 
user’s  choice.  Actually,  the  role  of  the  CWF  is  to  suppress  the  presence  of  terms  with  low 
order  derivatives,  which  are  responsible  for  destroying  the  strict  convexity  of  the  Laplace 
operator  Therefore,  because  of  strict  convexity,  rapid  global  convergence  of  a  number  of 
minimization  algorithms  to  the  unique  minimizer  of  J  is  guaranteed.  If  can  be  also  proven 
that  the  distance  between  the  minimizer  of  J  and  the  vector,  which  corresponds  to  the 
correct  solution  of  the  inverse  problem,  is  small,  if  the  data  are  not  ’’too  noisy”.  Preliminary 
results  in  this  direction  are  quite  promising. 
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